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On the basis of the Extensive Air Shower (EAS) data observed by the GAMMA experiment the en- 
ergy spectra and elemental composition of the primary cosmic rays have been derived in the 10 3 -i-10 5 
TeV energy range. Reconstruction of the primary energy spectra is carried out in the framework 
of the SIBYLL and QGSJET interaction models and the hypothesis of the power-law steepening 
primary energy spectra. The obtained energy spectra of primary H, He, O, Fe nuclei along with the 
SIBYLL interaction model agree with the corresponding extrapolations of known balloon and satel- 
lite data at the ~ 10 3 TeV energies. The energy spectra obtained from the QGSJET model, show 
predominant proton composition of cosmic rays in the knee region. The evident rigidity-dependent 
behavior of the primary energy spectra for both interaction models are displayed at the following 
rigidities: Er ~ 2400 -r 3000 TV (SIBYLL) and E R ~ 3400 ± 200 TV (QGSJET). 
Using parametric event-by-event method of the primary energy evaluation by measured 
N c h, N l± (E ti > 5GeV,R < 50m) and age (s) shower parameters, the all-particle energy spectra 
were obtained. 

All presented results are derived taking into account the detector response, reconstruction uncer- 
tainties of EAS parameters and fluctuation of EAS development. 

PACS numbers: 96.40.Pq, 96.40.De, 98.70.Sa 



I. INTRODUCTION 



The investigation of the energy spectra and elemental 
composition of primary cosmic rays in the knee region 
(10 3 -^ 10 5 TeV) remains to be one of the intriguing prob- 
lems of the modern high energy cosmic-ray physics. De- 
spite the fact that these investigations have been carried 
out for more than half a century, the data on the elemen- 
tal primary energy spectra at energies of E > 10 3 TeV 
need improvement. However, a bend of the all-particle 
energy spectra at around 3 • 10 3 TeV (called the "knee") 
at overall spectrum ~ E~ 27 until the knee and ~ E~ 31 
beyond the knee, may be considered as an avowed fact. 
Moreover, assuming that the supernova explosion is the 
main source of the cosmic rays, different theoretical mod- 
els of the high energy cosmic-ray origin and propagation 
through the Galaxy, predict the rigidity-dependent steep- 
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ening primary energy spectra in the knee region 0, 0, 0] . 
High statistical accuracies of the last EAS experiments 
[10 El already allowed us to infer that the rigidity- 
dependent steepening energy spectra of primary nuclei 
can approximately describe the observed EAS size spec- 
tra in the knee region in the framework of conventional 
interaction models. However, the accuracies of the ob- 
tained elemental primary energy spectra are still insuf- 
ficient due to both the uncertainty of interaction model 
and the accuracy of the solutions of the EAS inverse prob- 
lem. 

The Gamma facility (Fig. 1) was designed at the begin- 
ning of 90's in the framework of the ANI experiment 
and the preliminary results of EAS investigations pre- 
sented in H 0, 0, ^| . The main characteristic features 
of the GAMMA experiment are the mountain disposi- 
tion, symmetric location of the EAS detectors and under- 
ground muon scintillation carpet that detects EAS muon 
components at E^ > 5 GeV energies. 
Here, the description of GAMMA facility, EAS inverse 
approach determining the primary energy spectra using 
the observed EAS data, and main results of investigation 
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during 2002-2004 0, [n| are presented in comparison 
with the corresponding' MC-simulated data in the frame- 
work of the SIBYLL Hi and QGSJET [H interaction 
models. 



II. GAMMA EXPERIMENT 

The GAMMA installation is a ground based array of 33 
surface particle detection stations and 150 underground 
muon detectors located on the south side of Mount Ara- 
gats, Armenia. Elevation of the GAMMA facility is 3200 
m above sea level, which corresponds to 700 g/cm 2 of 
atmospheric depth. The diagrammatic layout is shown 
in Fig. 1. 

The surface stations of the EAS array are located on 5 
concentric circles of radii: 20, 28, 50, 70, 100 m and each 
station contains 3 square plastic scintillation detectors 
with the following dimensions: 1x1x0.05 m 3 . Each of the 
central 9 stations contains an additional (4-th) small scin- 
tillator with dimensions 0.3x0.3x0.05 m 3 (Fig. 1) for high 
particle density (3> 10 2 particles/m 2 ) measurements. 
The photomultiplier tube is positioned on the top of the 
aluminum casing covering the scintillator. One of the 
three station's detectors is examined by two photomulti- 
pliers, one of which is designed for fast-timing measure- 
ments. 

150 underground muon detectors (muon carpet) are 
compactly arranged in the underground hall under 2.3 
Kg/cm 2 of concrete and rock. The dimensions, casing 
and applied photomultipliers are the same as in the EAS 
surface detectors. 



A. Detector system and triggering 

The output voltage of each photomultiplier is con- 
verted into the pulse burst by logarithmic ADC and 
transmitted to the CAMAC array where the correspond- 
ing electronic counters produce a digital number (" code" ) 
of pulses in the burst. Four inner ("trigger") stations are 
monitored by a coincidence circuit. If each of at least 
two scintillators of each trigger station detects more than 
3 particles, the information from all detectors are then 
recorded along with the time between the master trigger 
pulse and the pulses from all fast-timing detectors. The 
given trigger condition provides EAS detection with the 
EAS size threshold N ch > (0.5 4- 1) • 10 5 at the location 
of the EAS core within the R <C 25 m circle. 
Before being placed on the scintillation casing, all pho- 
tomultipliers are tested by a test bench using a lumin- 
odiode method where the corresponding parameters of 
logarithmic ADC and the upper limits ((0.5 -f- 1) • 10 4 ) 
of the measurement ranges are determined. The num- 
ber of charged particles (m) passing through the i-th 
scintillator is computed using a logarithmic transforma- 
tion: hxrii = (C — Co)/d, where the scale parameter 
d ~ (9 -7- 10) ±0.35 is preliminarily determined by the test 



bench, C — (042 7 — 1) is an output digital code from the 
CAMAC array corresponding to the energy deposit of n 
charged particles into the scintillator, Co — (5 4 6) ±0.25 
is determined for each hour of run and is equal to the 
mode of the background single particle digital code spec- 
tra (Fig. 2). 

The time delay At, = tj — t\ of each j-th (j — 
2, . . . ,33) fast timing detector is estimated by the pair- 
delay method at the resolution time about 4 4 5 ns. 



B. Reconstruction of EAS parameters 

EAS zenith angle (8) is estimated on the basis of mea- 
sured shower front arrival times by 33 fast-timing sur- 
face detectors, applying the maximum likelihood method 
and flat-front approach |14L Il5j . The corresponding un- 
certainty are tested by MC simulations and is equal to: 
a(6) ~ 1.5°. 

The reconstruction of the EAS size (N c h), shower age 
(s) and core coordinates (xo,yo) are performed based 
on the NKG approximation of measured charged par- 
ticle densities ({ni},z = l,...,m) using the x 2 min- 
imization to estimate xo,yo and the maximum likeli- 
hood method to estimate the N c h taking into account 
the measurement errors. The logarithmic transforma- 
tion L(rij) = lnrii — (l/m)^hinj at ^ 0, allows to 
obtain the analytical solution for the EAS age parameter 
(s) at the x 2 minimization 0,0]. Unbiased (< 5%) 
estimations of N c h, s, xo,yo shower parameters are ob- 
tained at N ch > 5 • 10 5 , 0.3 < s < 1.6, and R < 25 m 
from the shower core to the center of the EAS array dis- 
tances. Corresponding accuracies are derived from MC 
simulations by the CORSIKA(EGS) [HI and are equal 
to: AN ch /N ch ~ 0.1, As ~ 0.05, Ax, Ay -0.5^1 m. 
The reconstruction of the total number of EAS muons 
(Nfj,) by the detected muon densities ({n^j},j = 
1, . . . , 150) from the underground muon hall is carried out 
by restricting the distance i? M < 50 m from the shower 
core (so called "truncated" EAS muon size [20j ) and 
Greisen approximation of the muon lateral distribution 
function: p^(r) = cN^(R < 50m) exp (— r/r )/(r/ro) ' 7 , 
where ro = 80 m, c = 1/2-7T J Q p(r)rdr. The trun- 
cated muon size N^(R < 50m) is estimated at known 
(from the EAS surface array) shower core coordinates 
in the underground muon hall. The unbiased estima- 
tions for muon size are obtained at > 10 3 using 
the maximum likelihood method and assuming Poisson 
fluctuations of detected muon numbers. The reconstruc- 
tion accuracies of the truncated muon size are equal to 
AN^/N^ ~ 0.2 ~ 0.35 at % ~ 10 5 ~ 10 3 respectively. 
Notice, that the detected muons in the underground hall 
are always accompanied by the electron-positron equi- 
librium spectrum which is produced when muons pass 
through the matter (2300 g/cm 2 ) over the scintillation 
carpet. Since this spectrum depends on the muon energy 
(~ hxE^), overestimations (~ 25 ~ 30%) of the recon- 
structed muon size have to depend on the primary energy 
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FIG. 1: Diagrammatic layout of the GAMMA facility 



and therefore on the EAS size. 



III. EAS SIMULATION 

A. Key assumptions 

All observed quantities (AF/Aq u ) in the high en- 
ergy EAS physics are obtained via convolutions of 
the energy spectra (Ha/cLE of primary nuclei (A 
II. He, ... at least up to Ni) with the differential spec- 
tra Wa(E, q u ) of the EAS parameters q u = N c h, N^, s at 
the observation level and EAS array response functions 
dK A (E,q u ,8)/dq u : 

af u i ^ r di A r r dn A 

(i) 

where the EAS parameter q u is a reconstructed value of 
the corresponding parameter q u on the observation level, 
dD = cos Odxdydft is an element of the multidimensional 
phase space (D) of the EAS detection taking into ac- 
count the EAS selection criteria and trigger conditions, 
W A (E, q u ,9) are the corresponding differential spectra 
of the EAS parameters (q u ) at the primary energy E, 
zenith angle of incidence 9 and a given kind of primary 
nucleus (A), C is a corresponding normalization factor. 
In t he g eneral case, Wa depends on the interaction model 

mm 

The multidimensional integral above is better to calcu- 
late by Monte-Carlo simulation, especially since the spec- 
tra Wa(E, q u , 9) can be computed more or less precisely 
only by 3-dimensional EAS simulations. 

B. Simulation scenario 

We have computed the shower spectra dWA(E, q u , 9), 
(q u = N c h, Np, s . . .) on the observation level of the 



GAMMA facility usingthe CORSIKA6031(NKG,EGS) 
EAS simulation code O with the QGSJET01 [l]| and 
SIBYLL2.1 ^3 interaction models for 4 groups (A = 
H, He, O, Fe) of primary nuclei at the power-law energy 
spectra (~ E~ 1 - 5 ) in the 5T0 2 -^5T0 5 TeV energy range. 
The spectral index (-1.5) was chosen to provide high sta- 
tistical accuracies of the simulated data beyond the knee. 
The EGS mode of the CORSIKA was used for computa- 
tions of the response functions of the GAMMA detectors 
taking into account the EAS gamma-quanta contribu- 
tions and the choice of the corresponding input parame- 
ters of the adequate NKG mode. 

Each EAS particle (7, e, fi, h) obtained from the COR- 
SIKA(EGS) on the observation level (not interrupting 
the CORSIKA routine) is passing through the steel cas- 
ing (1.5 mm) of detector station and then through the 
corresponding scintillator. The pair production and 
Compton scattering processes are additionally simulated 
in the case of the EAS 7-quanta passing through the steel 
casing and the scintillator. 

The resulting energy deposit in the scintillator is con- 
verted to the ADC code and inverse decoded into a num- 
ber of "detected" charge particles taking into account 
all uncertainties of the ADC parameters (Co , d) and the 
fluctuation of the light collected by the photomultiplier 
{01 ~ 0.25/ y/n). 

Using the simulation scenario above, 100 EAS events 
were simultaneously simulated by the CORSIKA rou- 
tine at the EGS and NKG modes for A = H, He, O, Fe 
primary nuclei at log-uniform energy spectra in the 
5 • 10 2 -i- 10 5 TeV energy range. The computations of 
the charged particle densities in the surface detectors at 
NKG mode of the CORSIKA were performed by applying 
the two-dimensional interpolations of the corresponding 
particle density matrix from the CORSIKA routine [13. 
The agreement (~ 5%) of the EGS and NKG simulated 
data was attained at the E e ~ 1 ± 1 MeV kinetic en- 
ergy threshold of the EAS electrons (positrons) at NKG 
mode (input parameter of CORSIKA code). However, 
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the energy threshold for the detection of a vertical sin- 
gle minimal ionizing background particle by scintillation 
counters is about 8-^9 MeV and differences obtained by 
the CORSIKA NKG prediction are completely explained 
by contribution of EAS 7-quanta. 

Thus, the EAS simulations by the CORSIKA with fast 
computation NKG mode at E e > 1 MeV is adequate to 
the EAS simulation by the EGS mode taking into ac- 
count the EAS 7-quanta and peculiarity of the GAMMA 
surface array. 

All EAS muons with energies of > 4 GeV on the 
GAMMA observation level have passed through the 2.3 
Kg/cm 2 of rock to the muon scintillation carpet (the un- 
derground muon hall) . Fluctuations of the muon ioniza- 
tion losses and electron (positron) accompaniment due to 
the muon bremsstrahlung, direct pair production, knock 
on and photo- nuclear interactions are taken into account. 
The transformation of the energy deposit to the number 
of detected muons is performed the same way as for the 
surface detectors. 

The EAS simulations were performed at 4.5- 10 4 primary 
H, 4.3 • 10 4 He, 2.4 • 10 4 O, 2.4 • 10 4 Fe nuclei using 
the CORSIKA NKG routine at the SIBYLL interaction 
model. Corresponding statistics at the QGSJET inter- 
action model were: 4.1 • 10 4 , 4.2 • 10 4 , 2 • 10 4 , 2 • 10 4 . 
The energy thresholds of primary nuclei were the same 
for both interaction models and were set 0.5, 0.7, 1, 1.2 
PeV respectively at 5 • 10 3 PeV upper energy limit. 



IV. MEASUREMENT ERRORS AND DENSITY 
SPECTRA 

The close disposition of k = 1,2,3 scintillators in each 
of the (z-th) detector station of the GAMMA surface ar- 
ray allows to auto-calibrate the measurement error by 
detected EAS data. The measured and simulated parti- 
cle density divergences (rife — p)/p versus average value 
p = (l/3)J2 n k at Ri > 10 m distances from shower 
core are shown in Fig. 2 (circle symbols). The obtained 
dependences are completely determined by Poisson fluc- 
tuations (at Ri 1 m ) and measurement errors. 

The agreement of the measured and simulated depen- 
dences allowed to extract the real measurement errors of 
the GAMMA detectors. In Fig. 2 the corresponding re- 
sults are shown (square symbols). 

The background single particle spectra (in the units of 
ADC code) detected by GAMMA surface scintillators 
for 78 sec operation time are shown in Fig. 3 (dotted 
lines). The background single particle spectra detected 
by underground muon scintillators have the same shape 
at about 10 times less intensities. 

These spectra are used for the operative (each hour) de- 
termination of ADC parameters (Co) during an exper- 
iment. The symbols and solid lines in Fig. 3 display 
the corresponding expected spectra obtained by MC- 
simulation taking into account the measurement errors 
(symbols) and without errors (line) respectively. The 
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FIG. 2: Particle density divergences (circle symbols) and 
measurement error of single detector (square symbols) ver- 
sus charged-particle density. 
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FIG. 3: Background single particle spectra of 15 surface detec- 
tors (dotted lines). The symbols (solid line) are the expected 
spectra taking into account (without) measurement errors. 



minimal primary energy in simulation of the background 
particle spectra was confined to the 7.6 GV primary par- 
ticle's geomagnetic rigidity. 

Because the effective primary energies responsible for the 
single particle spectra at observation level 700 g/cm 2 are 
about ~ 100 GeV and the energy range is studied by 
direct measurements in the balloon and satellite experi- 
ments, the primary energy spectra and elemental compo- 
sition at MC-simulation were taken from approximations 
plf. Notice, that the expected single particle spectra at 
these energies are practically the same for QGSJET and 
SIBYLL interaction models. 

Fig. 4a,b (symbols) display the charged particle density 
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spectra detected by the corresponding surface detectors 
(a) and underground muon detectors (b) at Ri < 50 
m and different EAS size thresholds: N ch > 5 • 10 5 , 
N ch > 10 7 (and additionally N ch > 2 • 10 6 for muon 
density spectra). 

The showers were selected at 9 < 30° and the shower 
core location in the R < 25 m range from center of the 
GAMMA facility (Fig. 1). The corresponding expected 
spectra (lines) at different interaction models are also 
shown in Fig. 4. The primary energy spectra and ele- 
mental composition at MC-simulations were taken from 
EAS inverse problem solution (see below). There is a 
good agreement of the expected and observed data for 
the surface array in the measurement range (about four 
orders of magnitude). However, the agreement of the 
detected muon density spectra with expected ones is at- 
tained only in the N c h < 10 7 range. The observed dis- 
crepancies for the muon density spectra at N c h > 10 7 
are unaccounted for the present and demand subsequent 
investigations. 

V. EAS DATA 

The main EAS data of the GAMMA experiment are 
shown in Fig. 5-10 (symbols). These results were ob- 
tained at the 6.19 • 10 7 sec operation time and following 
selection criteria: N ch > 5 • 10 5 , R < 25 m, 6 < 30°, 
0.3 < s < 1.6. All the lines and shaded areas in Fig. 5- 
10 correspond to the expected spectra according to the 
QGSJET and SIBYLL interaction models. 
The EAS size spectra (N^ h 5 ■ dF{6)/dN ch ) at 3 zenith 
angular intervals are shown in Fig. 5. The truncated 
muon size spectra in the same zenith angular intervals are 
shown in Fig. 6. These spectra normalized to the EAS 
intensity at N ch > 5-10 5 and 9 < 30°. The EAS size spec- 
tra at 9 < 30° and different thresholds of the truncated 
EAS muon size are shown in Fig. 7. The normalized EAS 
truncated muon size spectra at different EAS size thresh- 
olds are shown in Fig. 8. Fig. 9 displays the average EAS 
age parameter dependence on EAS size. The lines are 
the expected dependences according to QGSJET (dot- 
ted line) and SIBYLL (solid line) models. The obtained 
N^Nch) dependences and corresponding expected values 
at the primary Hydrogen, Iron and mixed compositions 
computed in the frame of the SIBYLL and QGSJET in- 
teraction models are plotted on Fig. 10. 



VI. EAS INVERSE PROBLEM AND PRIMARY 
ENERGY SPECTRA 

A. Combined approximations of EAS data 

Direct computations of the expected EAS spectra us- 
ing the integral expression (1) is possible only in the 
framework of a given interaction model and known pri- 
mary energy spectra. Moreover, the Gamma data shown 



in Fig. 4-10 may only formally compare with the same 
data obtained by other EAS experiments performed at 
both similar atmospheric depths [lj| |22, an d depths 
close to the sea level 0, 0, The correct compar- 

ison is possible only at known primary energy spectra 
and known interaction model because both transforma- 
tion of the detected EAS spectra to the spectra at a given 
observation level and the extrapolation of the obtained 
spectra to an another atmosphere depth in a general case 
are folded by the integral expressions similar to (1). 
In such case the more reliable way to interpret the ex- 
perimental data is to unfold the integral expression (1) 
at a given interaction model. As a criterion of the valid- 
ity of the solutions, the x 2 test of the detected and ex- 
pected data may be performed. The agreement between 
the obtained energy spectra at different primary nuclei 
and the corresponding extrapolations of known balloon 
and satellite data to the given measurement range will 
also validate the solutions. 

Evidently, the accuracies of the unfolding of expression 
(1) depend not only on number of measurement points 
(bins) and different measured spectra but also on the 
wealth of information about the primary energy spectra 
and the interaction model involved in the given measured 
EAS spectra. The amount of information contained in 
the expression (1) reveals itself via stability and uncer- 
tainties of the solutions. 

It is shown in 26], that the EAS size spectra and EAS 
truncated muon size spectra at three zenith angular in- 
tervals allow to reliably unfold expression (1) at a given 
interaction model for not more than 2 kinds of primary 
nuclei. The unreliability of solutions of (1) for 4 kinds of 
primary nuclei was shown in |27j . as well. 
Taking into account the above, we used the parameter- 
ization of the integral equation (1) similar to 0, l28| . 
The solutions for the primary energy spectra in (1) were 
sought based on the theoretically known power-law func- 
tion with the "knee" at the rigidity-dependent ener- 
gies Ek{A) = Er ■ Z and the same indices (—71) and 
(—72) before and after the knee respectively, for all kinds 
of primary nuclei (A): 




where 7 = 71 = 2.65 at E < Ek(A), 7 = 72 at 
E > Ek{A), Er is particle's rigidity and Z is a charge of 
A nucleus. 

Thus, the integral equation (1) is transformed into a 
parametric equation with unknown spectral parameters: 
<£(A), Ek(A), 72, which are determined by minimization 
of x 2 function: 

2 _ 1 V - ^ (Cu,v — S,u,v) 2 /„n 

where U is the number of examined functions 
Q u ^ v = AF U / Aq u ^ v (Fig. 5-10, symbols) obtained 
from the experiment with statistical accuracies cr(£ UtV ) 
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FIG. 4: Detected (symbols) and expected (lines) particle density spectra of surface scintillators (left panel, a) and underground 
muon scintillators (right panel, b). 




FIG. 5: EAS size spectra at 3 zenith angular intervals 
(symbols) and corresponding expected spectra according to 
the SIBYLL (shaded area) and QGSJET interaction models 
(lines) . 



FIG. 6: Normalized EAS truncated muon size spectra at 3 
zenith angular intervals (symbols). The lines correspond to 
the expected spectra at the SIBYLL (solid) and QGSJET 
(dashed) interaction models. 



at v — 1,...,V U measured points (bins), and £ UjV and 
are the corresponding expected values of the 
examined data set. 

Using the aforementioned formalism and U — 6 2- 
dimensional examined functions from Fig. 5-8 (symbols) 
and 1-dimcnsional functions from Fig. 9,10, the unknown 
spectral parameters <&(A), Ek(A), 72 were derived by the 
minimization of x 2 function (3) at 71 = 2.65 and the de- 
gree of freedom V u ~ 350. 

The values of spectral parameters (2) obtained by the so- 
lution of the parameterized equation (1) are presented in 



Table 1 at the QGSJET and SIBYLL interaction mod- 
els. The derived primary energy spectra for p, He, O, Fe 
nuclei are shown in Fig. 11 (shaded areas) in comparison 
with the KASCADE data (symbols) from p5| . 
The expected spectra conforming the examined data set 
according to the solutions above are shown in Fig. 5-10 
(lines and shaded area) for the QGSJET and SIBYLL 
interaction models. It is necessary to note, that the ob- 
tained results in the framework of the SIBYLL interac- 
tion model are more consistent and slightly dependent on 
a number of examined functions. 
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FIG. 7: EAS size spectra (symbols) at different truncated 
muon size thresholds and 8 < 30°. The lines and shaded 
areas are the expected spectra according to the QGSJET and 
SIBYLL interaction models. 



FIG. 9: Dependence of average EAS age parameter on EAS 
size. 
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FIG. 8: Normalized EAS truncated muon size spectra at dif- 
ferent EAS size thresholds (symbols). The lines and shaded 
areas correspond to the expected spectra at the QGSJET and 
SIBYLL interaction models respectively. 



B. 4-Dimensional approach 



The combination of 1,2-dimensional approximations of 
EAS data above does not take into account all the in- 
formation about primary energy spectra folded in the 
detected EAS data. In general, the EAS inverse problem 
can be formulated in the multidimensional space of EAS 
parameters. In case of the 4-parametric (N c h, N^, s,0) 
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FIG. 10: Average EAS truncated muon size N^(E > 5 GeV, 
R < 50m) versus EAS size (N c h). The dashed lines and 
shaded area are the expected dependences at the SIBYLL 
interaction models and pure p, Fe and mixed compositions 
respectively. Dash-dotted and solid lines are the same depen- 
dences at the QGSJET models. 



analysis, the expression (1) is written as: 



AF 



AN ch ANuASAil 



xTZ(9)dEdDdN ch dN f ,ds , (4) 



where Ga{E,6) = d 3 W A {E,6)/dN ch dN tl ds, are the 
multidimensional differential EAS spectra at given 
A, E, 9 parameters of the primary nucleus, TZ(6) = 
d 3 R(0) /dNchdN^ds are the error functions of the ex- 
periment. The parameters with a tilde symbol are the 
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FIG. 11: Energy spectra and abundance of the primary nuclei (shaded areas) at the SIBYLL (left panel) and QGSJET (right 
panel) interaction models. The symbols are the KASCADE data from |25l| . 



TABLE I: Parameters of primary energy spectra (2) at 1,2- 
dimensional analysis of EAS data. Scale factors $a and par- 
ticle's rigidity Er have units of (m 2 ■ sec ■ ster ■ TeV) -1 and 
(TV) respectively. 



Parameters 



SIBYLL 



QGSJET 



<&He 

Er 
"71 

;r 



0.081 ±0.004 
0.072 ± 0.008 
0.028 ± 0.008 
0.028 ± 0.003 
2560 ± 200 

2.65 
3.21 ± 0.04 
2.5 



0.164 ±0.004 
0.005 ± 0.008 
0.005 ± 0.006 
0.018 ±0.003 
3400 ± 150 

2.65 
3.10 ±0.03 
2.6 



"Parameter was fixed. 



reconstructed values of corresponding EAS parameters. 
Evidently, the amount of information about primary en- 
ergy spectra contained in the detected multidimensional 
spectrum AF is always greater than the cumulative 
amount of information contained in the 1,2-dimcnsional 
spectra AF U /Aq u , q u = N c h, N^, s of the expression (1). 
The difference is determined by the inter-correlations of 
EAS parameters that are taken into account in the ex- 
pression (4). 

On the basis of the EAS data set of the GAMMA exper- 
iment, the simulated EAS database (section III) and pa- 
rameterization (2), the equations (4) were resolved by the 
X 2 -minimization method. The computations were per- 
formed at the following bin dimensions: A In N c h = 0.15, 
A In = 0.25, Asec6> = 0.05 and As = 0.15 on the 
left and right hand side of s* = 0.85 and As = 0.3 in 
other cases. The total number of the degree of freedom 
at 4-dimensional ^-minimization was equal to 1560. 
The values of spectral parameters (2) obtained by the so- 



TABLE II: Parameters of primary energy spectra (2) at 4-D 
analysis of EAS data. Scale factors 3>a and particle's rigidity 
Er have units of (m 2 • sec- ster -TeV) -1 and (TV) respectively. 

Parameters SIBYLL QGSJET 



<E>o 

Er 

72 



0.089 ± 0.003 
0.053 ± 0.005 
0.049 ± 0.004 
0.029 ± 0.003 
3000 ± 300 
3.16 ±0.08 
1.2 



0.14 ±0.004 
0.034 ± 0.004 
0.016 ±0.003 
0.015 ±0.002 
3300 ± 200 
3.10 ±0.03 
1.1 



lution of the parameterized equation (4) are presented in 
Table 2 at the QGSJET and SIBYLL interaction mod- 
els. As it is seen from Fig. 11 and Tables 1,2, the de- 
rived expected primary energy spectra significantly de- 
pend on interaction model. The expected abundance of 
primary nuclei at energy E ~ 10 3 TeV in the framework 
of SIBYLL model agrees well with corresponding extrap- 
olations of the balloon and satellite data [2l|, whereas 
the predictions according to the QGSJET model point 
out to a predominantly proton primary composition in 
the 10 3 ± 10 5 TeV energy range. 



VII. EVENT-BY-EVENT ANALYSIS 

The mountain location of the GAMMA experiment 
and the agreements of observed and simulated data in 
the measurement range 5 • 10 5 < N c h < 10 7 (Fig. 4-10) 
allowed, apart from above, to obtain the all-particle en- 
ergy spectra with high reliability . The method is based 
on an event-by-event evaluation of primary energy us- 
ing the reconstructed parameters N c h,Nfi, s, 6 of detected 
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EAS. Such possibilities have been studying for a long 
time in different papers 0, |3(| and the main dif- 
ficulty was to obtain an unbiased energy estimation at 
an existent abundance of the primary nuclei taking into 
account the fluctuations of shower development and de- 
tector response. 

Using the simulated database, J — 1.5 ■ fO 4 EAS events 
were taken for each of k = 1, . . . , 4 kinds (H, He, O, Fe) 
of primary nuclei and each interaction model (S1BYLL, 
QGSJET). The reconstructed N c h,Na,s shower param- 
eters, known zenith angle 9 and primary energy E were 
used at minimization 



, a p ) = 



1 4 - 



(fn.£a,fc,j _ ]nEo,k,jY 



(5) 



fc=i j=i 



where E\ — f(a\, . . . , a p \N c h, N^, s, 9) is the investigated 
parametric function with a±, . . . , a p parameters, cje is ex- 
pected accuracy of the E\ evaluated energy. The best 
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FIG. 12: Accuracy (RMSD) of the primary energy evaluations 
at different number of approximation parameters. 



estimations were achieved at 7-parametric (p = 7) fit: 

:, (6) 



\nEi=aiX-\ h a 3 + 04c + a 5 e + .. 

c (x - a 7 y) • 



where x = lnN c h, y — In N^(R < 50m), c = cos 9 and en- 
ergy Ei has units of GeV. The values of the oi, . . . , a 7 pa- 
rameters for both interaction models and the correspond- 
ing x 2 obtained from (5) at erg = 0.15 are displayed in 
Table 3. The root mean square deviations of the energy 
estimation by 7-parametric fit (4) in the framework of 
the SIBYLL model is shown in Fig. 12. The correspond- 
ing results at three (only x, s variables) and 4-parametric 
(x, s, cos 9) fit are shown in Fig. 12 as well. 
The obtained error distributions estimating primary en- 
ergy by 7-parametric approximation (6), are shown in 
Fig. 13 for H, He, O, Fe nuclei. The red line corresponds 



TABLE III: Approximation parameters oi, . . . , 07 of primary 
energy evaluation (6) and corresponding \ 2 obtained from (6) 
at the SIBYLL and QGSJET interaction models. 



Model 


ai 


a 2 


as 


"4 


a-5 


a 6 


a 7 


x 2 


SIBYLL 


1.03 


3.98 


-4.3 


2.01 


-1.2 


11.8 


0.94 


0.85 


QGSJET 


1.03 


4.38 


-4.6 


2.35 


-1.3 


11.5 


0.96 


0.94 
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FIG. 13: Distribution of errors of the primary energy estima- 
tion by event-by-event 7-parametric fit at different primary 
nuclei. 



to the Gaussian distribution at the same parameters as 
cumulative distribution (black solid line). 
The all-particle energy spectrum derived on the basis of 
the GAMMA 2002-2004 EAS data set and fit (6), at 
the QGSJET (filled red square symbols) and SIBYLL 
(filled blue circle symbols) interaction models are shown 
in Fig. 14. 

Notice, that the energy spectrum obtained by event- 
by-event method claims additional corrections, because 
the errors <te = a(AE/E) and power-law energy spec- 
tra (~ E^ 1 ) lead to an overestimation of the spectrum 
r\ = exp (((7 — 1)cte) 2 /2) times. Moreover, the inevitable 
biases of energy estimations e(A) =< E\/Eq > depend 
on primary nuclei and shift the corresponding energy 
spectra /3(A) = e 1 ^ 1 times. The spectral shift due to 
(3(A) ^ 1 impossible to take into account without infor- 
mation about abundance of primary nuclei. 
The observed biases of 7-parametric fit (6) are distributed 
from e(p) ~ 1.02% up to e(Fe) ~ 0.96% (Fig. 13) and 
here are neglected. In the results shown in Fig. 14, the 
corrections of ri(E) are taken into account using the ex- 
pected accuracies from Fig. 12. 

The solid (red) and dashed (blue) lines in Fig. 14 rep- 
resent the all-particle primary energy spectra obtained 
on the basis of GAMMA data set by the solution of 
parametrized EAS inverse problem in the framework of 
the SIBYLL and QGSJET models respectively. The 
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15 o.9 - GAMMA data 

> 0.8 - SIBILL : • 7-fit, unfolded 

0.7 - QGSJET: "7-fit, ---unfolded 




Primary energy E (TeV) 

FIG. 14: All particle primary energy spectra obtained by 
event-by-event 7-parametric analysis (filled symbols) and 
EAS inverse problem solutions (solid and dashed lines) on 
the basis of GAMMA 2002-2004 database. 

event-by-event analysis of the GAMMA data at the 
QGSJET interaction model using a-parametric method 
[lfj| also shown in Fig. 14 (asterisk symbols). The dotted 
line in Fig. 14 represents the parametrized solutions of 
the EAS inverse problem for the KASCADE EAS data 
at rigidity-dependent steepening primary energy spectra 
[H. The results of KASCADE02 in Fig. 14 obtained 
by the non-parametric event-by-event analysis was taken 
from review H3- The KASCADE01,05 data obtained by 
the iterative method [33| of unfolding of the EAS inverse 
problem were taken from [23, respectively. 

VIII. CONCLUSION 

The self-consistency of results (Fig. 2-10) obtained by 
GAMMA experiment at least up to N c h — 10 7 and cor- 
responding predictions in the framework of hypothesis of 
the rigidity-dependent steepening primary energy spec- 
tra and validity of the SIBYLL or QGSJET interaction 
models point towards: 

• The anomalous behavior of the EAS muon spectra 



(overestimation, Fig. 4b, 8, 10) and EAS age param- 
eter at EAS size N c h > 10 7 . The same behavior of 
the EAS age parameter had been observed also in 

mm 

• The obtained abundances and energy spectra of 
primary p, He, O, Fe nuclei depend on interaction 
models. The SIBYLL interaction model is more 
preferable in terms of the extrapolation of the de- 
rived expected primary spectra (Fig. 11) to the en- 
ergy range of the direct measurements. 

• The rigidity-dependent steepening energy spectra 
of primary nuclei describe the EAS data of the 
GAMMA experiment at least up to N ch ~ 10 7 
with average accuracy < 10% at particle's mag- 
netic rigidity E R ~ 2400-^3000 TV (SIBYLL) and 
E R ~ 3300 ± 200 TV (QGSJET). 

• The 4-dimensional approach at the EAS inverse 
problem solution is more preferable in terms of the 
stability and accuracies of solutions. 

• The obtained all-particle energy spectra slightly 
depend on interaction model and are practically 
the same at both the event-by-event reconstruction 
method and the EAS inverse approach. 

The obtained energy spectra of primary nuclei {A = 
p, He, O, Fe) in the energy range 10 6 < Ea < 5 • 10 7 
GeV disagree (see Fig. 11) with the same KASCADE 
data obtained by the iterative method [3j|. However, the 
discrepancies of all-particle energy spectra (see Fig. 14) 
obtained by the GAMMA and KASCADE experiments 
are sufficiently small (~ 20%). 
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